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Abstract — Due to its reduced communication overhead and 
robustness to failures, distributed energy management is of 
paramount importance in smart grids, especially in microgrids, 
which feature distributed generation (DG) and distributed stor- 
age (DS). Distributed economic dispatch for a microgrid with high 
renewable energy penetration and demand-side management 
operating in the grid-connected mode is considered in this paper. 
To address the intrinsically stochastic availability of renewable 
energy sources (RES), a novel power scheduling approach is 
introduced. The approach involves the actual renewable energy as 
well as the energy traded with the main grid, so that the supply- 
demand balance is maintained. The optimal scheduling strategy 
minimizes the microgrid net cost, which includes DG and DS 
costs, utility of dispatchable loads, and worst-case transaction 
cost stemming from the uncertainty in RES. Leveraging the dual 
decomposition, the optimization problem formulated is solved 
in a distributed fashion by the local controllers of DG, DS, and 
dispatchable loads. Numerical results are reported to corroborate 
the effectiveness of the novel approach. 

Index Terms — Energy management, economic dispatch, de- 
mand side management, renewable energy, distributed energy re- 
sources, microgrids, robust optimization, distributed algorithms. 



I. Introduction 

Microgrids are power systems comprising distributed energy 
resources (DERs) and electricity end-users, possibly with con- 
trollable elastic loads, all deployed across a limited geographic 
area [TJ. Depending on their origin, DERs can come either 
from distributed generation (DG) or from distributed storage 
(DS). DG refers to small-scale power generators such as diesel 
generators, fuel cells, and renewable energy sources (RES), 
as in wind or photovoltaic (PV) generation. DS paradigms 
include batteries, flywheels, and pumped storage. Specifically, 
DG brings power closer to the point it is consumed, thereby 
incurring fewer thermal losses and bypassing limitations im- 
posed by a congested transmission network. Moreover, the in- 
creasing tendency towards high penetration of RES stems from 
their environment-friendly and price-competitive advantages 
over conventional generation. Typical microgrid loads include 
critical non-dispatchable types and also elastic controllable 
ones. 

Microgrids also entail distribution networks with residential 
or commercial end-users, in rural or urban areas. A typical 
configuration is depicted in Fig.[TJ see also The microgrid 
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Fig. 1 . Distributed control and computation architecture of a microgrid sys- 
tem. The microgrid energy manager (MGEM) coordinates the local controllers 
(LCs) of DERs and dispatchable loads. 

energy manager (MGEM) coordinates the DERs and the 
controllable loads. Each of the DERs and loads has a local con- 
troller (LC), which coordinates with the MGEM the scheduling 
of resources through the communications infrastructure. The 
microgrid can operate in two modes: either connected or 
disconnected from the grid, the latter referred to as island 
mode. 

In this context, the present paper deals with optimal energy 
management for both supply and demand of a microgrid 
incorporating RES. The microgrid is connected to the main 
grid, while energy can be sold to or purchased from the 
main grid. Decentralized algorithms are developed, which are 
robust to the uncertainty of the available RES. In addition to 
computationally efficient, distributed power scheduling must 
be resilient to communication outages or attacks. Furthermore, 
DER scheduling in microgrids must account for the variable 
and nondispatchable nature of the RES. 

Optimal energy management for microgrids including eco- 
nomic dispatch (ED), unit commitment (UC), and demand- 
side management (DSM) is addressed in J2), but without 
pursuing a robust formulation against RES uncertainty. Mixed 
integer programming problems are formulated for microgrid 
scheduling and DSM in (3), O- Based on deterministic RES 
models (e.g., those relating wind power with wind speed), 
ED problems are investigated in and (6|. In all the afore- 
mentioned works however, robust formulations accounting for 
the RES randomness are not pursued. Lyapunov stochastic 
optimization has been adopted to maximize the long-term 
profit of a RES facility in Q. Stochastic programming is 
also used to cope with the variability of RES. Single-period 
chance-constrained ED problems for RES have been studied 
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in JSJ, yielding probabilistic guarantees that the load will be 
served. Considering the uncertainties of demand profiles and 
PV power generation, a stochastic program is formulated to 
minimize the overall cost of electricity and natural gas for 
a building in j9)- Without DSM, robust scheduling problems 
with penalty-based costs for uncertain supply and demand 
have been investigated in IfTOl . Recent works explore energy 
scheduling with DSM and RES using only centralized algo- 
rithms ifTTl . Ifl2l . An energy source control and DS planning 
problem for a small microgrid system is formulated and 
solved using model predictive control in 1131 . Using energy 
storage to mitigate fluctuation in generation due to time- 
varying RES, an optimal power flow problem is formulated 
as a finite-horizon optimal control problem in 1141 . Game- 
theoretic approaches for microgrids using (non-)cooperative 
and dynamic potential games are introduced in lfl~5l and lfl6l . 
respectively. Distributed algorithms are developed in |fl7l . but 
they only coordinate DERs to supply a given load without 
considering the stochastic nature of RES. Similarly, robust 
energy management formulations have not been considered 
in the mentioned works lTT3l - lfl6l . 

This paper considers energy management of DERs and 
dispatchable loads with the goal of minimizing the microgrid 
net cost. The objective accounts for conventional DG cost, 
utility of elastic loads, penalized cost of DS, and a worst- 
case transaction cost. The latter stems from the ability of the 
microgrid to sell excess energy to the main grid, or to import 
energy in case of shortage. A robust formulation accounting 
for the worst-case amount of harvested RES is developed. A 
novel model is introduced in order to maintain the supply- 
demand balance arising from the intermittent RES. Moreover, 
a transaction-price-based condition is established to ensure 
convexity of the overall problem. The separable structure and 
strong duality of the resultant constrained problem are lever- 
aged to develop a low-overhead distributed algorithm based 
on dual decomposition. The distributed implementation relies 
upon message exchanges between the MGEM and LCs. For 
faster convergence, the proximal bundle method is employed 
for the non-smooth subproblem handled by the LC of RES. 

The rest of the paper is organized as follows. Section HI] 
formulates the robust energy management problem, and Sec- 
tion [HI] develops its distributed solver. Numerical results are 
reported in Section[lV] while conclusions and research outlook 
directions are provided in Section [V] 

Notation. Boldface lower case letters represent vectors; M" 
and R stand for spaces of n X 1 vectors and real numbers, 
respectively; M™ is the rt-dimensional non-negative orthant; 
x' transpose, and ||x|| the Euclidean norm of x. 



II. Robust Energy Management Formulation 

Consider a microgrid comprising M conventional (fossil 
fuel) generators, / RES facilities, and J DS devices (see also 
Fig. ID- The scheduling horizon is T := {1,2,..., T} (e.g., 
one-day ahead). The particulars of the optimal scheduling 
problem are explained in the next subsections. 



A. Load Demand Model 

Loads are classified in two categories. The first comprises 
inelastic loads, whose power demand should be satisfied at 
all times. Examples are power requirements of hospitals or 
illumination demand from residential areas. The total fixed 
power demand from these critical loads at slot t is denoted by 
L\ 

The second category consists of elastic loads, which are 
dispatchable, in the sense that their power consumption is 
adjustable, and can be scheduled. These loads can be further 
divided in two subclasses, denoted respectively by /Ci and /C2, 
each having the following characteristics: 

i) Class K,\ contains loads with power consumption £ 
[Pd™,Pd™], where n G N := {1,...,N}, and t e 
T- Higher power consumption yields higher utility for 
the end user. The utility function of the nth dispatchable 
load, U B {Pjj ), is selected to be strictly increasing and 
concave, with typical choices being piecewise linear or 
smooth quadratic. An example from this class is an A/C 
unit. 

ii) Class /C2 includes loads indexed by q 6 Q := {1, . . . , Q} 
with power consumption limits PJ? m and PJ? ax , and 
prescribed total energy requirements E q which have to 
be achieved from the start time S q to termination time 
T q . This type of loads can be the plug-in hybrid electric 
vehicles (PHEVs). Power demand variables {P E }J =1 

therefore are constrained as Y^t=s P e = ^1 P e e 
[P Eq ,P Eq '}, t€T, while P Eq =P Eq = for 
t ^ {S q , . . . ,T q }. Higher power consumption in earlier 
slots as opposed to later slots may be desirable for a cer- 
tain load, so that the associated task finishes earlier. This 
behavior can be encouraged by adopting for the qth load 
an appropriately designed time-varying concave utility 
function U Eq {P Eq ). An example is U Eq {P Eq ) := n q P Eq , 
with weights {71-*} decreasing in t from slots S q to T q . 
Naturally, U E (P E ) = can be selected if the consumer 
is indifferent to how power is consumed across slots. 

B. Distributed Storage Model 

Let Bj denote the state of charge (SoC) of the jth battery 
(i.e., the amount of energy storage) at the beginning of the 
slot t, with initial available energy B® while 5™ ax denotes 
the battery capacity, so that < Pj < B" lax , j e J := 
{1, . . . , J}. Let P B . be the power delivered to (drawn from) 
the jth storage device at slot t, which amounts to charging 
(Pg. > 0) or discharging (P| < 0) of the battery. Clearly, 
the SoC obeys the dynamic equation 

P* +1 =B) + P t B] , je J, teT. (l) 

Variables P B are constrained in the following ways: 
i) The amount of (dis)charging is bounded, that is 

-r,fB] < P B] < jjf {Bf " - B)) (3) 

with bounds Pg. in < and P^ ax > 0, while rjf and rff s 
are charging and discharging efficiencies; and 
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ii) Final SoC is also bounded, that is Bj > Bf m . 
The storage cost Pj(P*) can be also taken into account. 
Considering the battery reserve, a penalty can be postulated 
proportional to the deviation from the battery capacity as 
Hj(Bj) := n7*(P™ ax — Pj), meaning it is desirable to maintain 
a close-to-full battery level at the end of each period fl4l . 

C. Worst-case Transaction Cost 

Let W\ denote the actual renewable energy harvested by the 
ith RES facility at time slot t, and also let w collect all Wf, 
i.e., w := [Wl, W?, ...,W},..., Wj\. To capture the 
intrinsically stochastic and time-varying availability of RES, 
it is postulated that w is unknown, but lies in a polyhedral 
uncertainty set W. The following are two practical examples 
of uncertainty sets. 

i) The first example postulates a separate uncertainty set Wi 
for each RES facility in the form 



< e w i < w > 

teTi, 



max 



(4) 



where W\ (Wi) denotes a lower (upper) bound on W\\ 
time horizon T is partitioned into consecutive but non- 
overlapping sub-horizons %, s for i = I,..., I, s = 
1, 2, . . . , S; the total renewable energy for the ith RES 
facility over the sth sub-horizon is assumed bounded by 
W™ n and Wif*. In this example, uncertainty set W 
takes the form of Cartesian product 



W = Wi x . . . x Wi. 



(5) 



ii) The second example assumes a joint uncertainty model 
across all the RES facilities as 



from/to the spot market. Let be an auxiliary variable 
denoting the net power delivered to the microgrid from the 
renewable energy sources and the distributed storage in order 
to maintain the supply-demand balance at slot t. The shortage 

energy per slot t is given by 



while the surplus energy is 



where [a] + := max{a,0}, and [a]~ := max{— a, 0}. 

The amount of shortage energy is bought with known 
purchase price a*, while the surplus energy is sold to the 
main grid with known selling price /?*. The worst-case net 
transaction cost is thus given by 



(7) 



where {P^} collects P*. for t = 1,2, ...,T and {P*, } 
collects Pjj. for j = l,2,...,J, t = 1, 2, . . . , T. 

Remark 1. (Stochastic model). Instead of the worst-case 
robust model advocated here, a probabilistic approach is also 
possible. Specifically, suppose that W\ is a function of a 
random quantity v\. In the case of wind generation, v\ is the 
wind velocity, for which different models are available, and 
the mapping Wl(v\) is known |fl9l . Due to power shortage, 
the expected purchase cost becomes G({P^}, {P| .}) := 

ELi ^APk - Eli wf(vl) + zU where v * ; = 

[v 1 , . . . , v T ]. Likewise, the expected revenue can be obtained 
by averaging [P^ - Y! 1= i wf(vf) + £)/ =1 P%.]- over t. This 
stochastic model is useful when statistical knowledge for the 
RES is available, as when v stands for the wind velocity, 
which typically satisfies the Weibull distribution |fl9l . 



W := < w\W; <Wf < W, 



>T=\JT. 



(6) 



s=l 



where W_l (Wj) denotes a lower (upper) bound on W\\ 
time horizon T is partitioned into consecutive but non- 
overlapping sub-horizons T s for s = 1, 2, . . . , S; the total 
renewable energy harvested by all the RES facilities over 
the sth sub-horizon is bounded by I'F" 1111 and W r s max ; see 
also 072|. 

The previous two RES uncertainty models are quite 
general and can take into account different geographical 
and meteorological factors. The only information required 
is these deterministic lower and upper bounds, namely 
W$,W\,W™s n , ^m"' W™ in ,W™ a *, which can be deter- 
mined via inference schemes based on historical data [ 181- 



Supposing the microgrid operates in a grid-connected mode, 
a transaction mechanism between the microgrid and the main 
grid is present, whereby the microgrid can buy/sell energy 



D. Microgrid Energy Management Problem 

Apart from RES, microgrids typically entail also conven- 
tional DG. Let Pq be the power produced by the mth 
conventional generator, where m E M. := {1, ...,M} and 
(eT. The cost of the mth conventional generator is given by 
a strictly increasing and convex function C^Pq ). Typically, 
the chosen C* n (Pg ) is either piecewise linear or smooth 
quadratic. 

The energy management problem amounts to minimizing 
the microgrid social net cost; that is, the cost of conventional 
generation, storage, as well as the worst-case transaction cost 
(due to the volatility of RES) minus the utility of dispatchable 
loads: 



T / M 



N 



( pi > (P ™» t E E - E u k(phj 

\^G m >^D n < t=1 \ m =i „=1 



Q J \ 

E u k( p k) + E h K b j) + G (i pt R}> { p k }) 

9=1 J=l / 

(8a) 
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subject to: 



1 G„ 
pt 



G m 
M 



<n 
_ pt- 

G, 
- P* 



< Pg ax , m g M, teT 



1 < P. 



me M, teT 



Hn,up ) 

< -Rm.down, m e M, teT 



max 

G m 



pt 

G„ 



m—1 

pmin ^ pt 



^ pmax 



> SR', t e T 
neJV, teT 



/ n max.( _ , _ t~ 

<P En , qeQ, teT 



Y,P% q =E q , qeQ, 



t=S q 

o < b) < p; nax , 



< pi 



<s pmax 

^ "b 



dis pt < pt 



jej, teT 

, j ej, teT 
< v f(B™ 



(8b) 
(8c) 
(8d) 

(8e) 

(8f) 
(8g) 

(8h) 

(8i) 
(8j) 



B t +1 



Pi, jej, teT 



> Bf m , 



teT 



P 



<Pi<Pi 



m— 1 



pt 

^R 



ax , teT 

N 

Lt + j^pi 

n=l 



B)), jej, teT 
(8k) 

(81) 

(8m) 

(8n) 



q=l 



t e T. 



Constraints d8bl-(|8e]i stand for the minimum/maximum 
power output, ramping up/down limits, and spinning reserves, 
respectively, which capture the typical physical requirements 
of a power generation system. Constraints (f8fb and ([8nl 
correspond to the minimum/maximum power of the flexible 
load demand and committed renewable energy. Constraint ( f8oT > 
is the power supply-demand balance equation ensuring the 
total demand is satisfied by the power generation at any time. 

Note that constraints (f8bb— (f8ob are linear, while C^(-), 
— Ujj (•), -ffj(')' an d —UE q (-) are convex (possibly non- 
differentiable or non-strictly convex) functions. Consequently, 
the convexity of (PI) depends on that of G({P R }, {P B .}), 
which is established in the following proposition. 

Proposition 1. If the selling price (3* does not exceed the pur- 
chase price a* for any teT, then the worst-case transaction 
cost G({P B }, {P B .}) is convex in {P R } and {P Bj }- 

Proof: Using that [a} + + [a]~ = |a|, and [a} + — [a]~ = a, 
G({P R }, {P Bj \) can be re-written as 



G({j* }l {iV)=«u«$:[* 



P t 

r R 



i=l 



3=1 



+ l t (P, 




(9) 



with 5 t := (a* - /3')/2, and 7 * := (a 1 + (3 t )/2. Since 
the absolute value function is convex, and the operations 
of nonnegative weighted summation and pointwise maximum 
(over an infinite set) preserve convexity |20l Sec. 3.2], the 
claim follows readily. ■ 



An immediate corollary of Proposition Q] is that the energy 
management problem (PI) is convex if /3* < a* for all t. 
The next section focuses on this case, and designs an efficient 
decentralized solver for (PI). 

III. Distributed Algorithm 

In order to facilitate a distributed algorithm for (PI), a vari- 
able transformation is useful. Specifically, upon introducing 
Pr : = Pr. + E/=i p b^ ( pi ) can be re-written as 



N 



(P2) min£ [J2 CUPgJ- £ K^D 



t=l \m=l 



q=l j=l ) 

subject to: ([8bJ - d8o]i 

J 

pt R = PR + J2 Pt B^ teT ( 10b > 

3=1 

where x collects all the primal variables 
{P^P^ptPzB^PlP^y, {P R } collects P* 
for t = 1,...,T; and 



(8o) G({PA}) 



max 



+7 (Pk-J2 w ') 



(id 



The following proposition extends the result of Proposition [T] 
to the transformed problem, and asserts its strong duality. 

Proposition 2. If (P2) is feasible, and the selling price /?' 
does not exceed the purchase price a* for any teT, then 
there is no duality gap. 

Proof: Due to the strong duality theorem for the optimiza- 
tion problems with linear constraints (cf. [21, Prop. 5.2.1]), 
it suffices to show that the cost function is convex over the 
entire space and its optimal value is finite. First, using the 
same argument, convexity of G({P R }) in {P R } is immediate 
under the transaction price condition. The finiteness of the 
optimal value is guaranteed by the fact that the continuous 
convex cost ( 11 Oat is minimized over a nonempty compact set 
specified by (IFrJb— (l8ob. and ( llObb . ■ 
The strong duality asserted by Proposition |2] motivates the 
use of Lagrangian relaxation techniques in order to solve 
the scheduling problem. Moreover, problem (P2) is clearly 
separable, meaning that its cost and constraints are sums 
of terms, with each term dependent on different variables, 
namely {P^J, {P*J, {P|J, {Pj}, {P R } and {P R }. The 
features of strong duality and separability then imply that 
Lagrangian relaxation and dual decomposition are applicable 
to yield a decentralized algorithm; see also related techniques 
in the power systems ll22l and communication networks l23l . 
Coordinated by dual variables, the dual approach decomposes 
the original problem into several separate subproblems that 
can be solved by the LCs in parallel. The development of the 
distributed algorithm is undertaken next. 
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A. Dual Decomposition 



Constraints dSel l, ( f8ob , and dlObt couple variables across 
generators, loads, and the RES. Let z collect dual variables 
{/!*}, {A 4 }, and jV*}, which denote the corresponding La- 
grange multipliers. Keeping the remaining constraints implicit, 
the partial Lagrangian is given by 



M 



N 



t=l \m=l n=l 

Q J \ 

9=1 i=i / 

+e{^ ( SRt -f2( p ST- p Gj) 

t=l { \ m=\ ) 

N Q \ 

n=l q=l ) 



Then, the dual function can be written as 



(12) 



P(z) =min£(x, z) 

X 

s.t. <[8b]i - ([Hdll, © - <Hn) 



and the dual problem is given by 



max {A<}><}) 

s.t. /i* > o, A*,i/ el, t e T. 



(13a) 
(13b) 



The subgradient method will be employed to obtain the 
optimal multipliers and power schedules. The iterative process 
is described next, followed by its distributed implementation. 



1) Subgradient Iterations: The subgradient method 
amounts to running the recursions ll24l Sec. 6.3] 

(14a) 
(14b) 
(14c) 



fi^k + l) = [^{k)+ag^{k)] + 
X t {k + l) = X t (k)+ag xt (k) 
v l {k + 1) = v\k) + ag u t(k) 



where k is the iteration index; a > is a constant stepsize; 
while g^t (fc), g\± {k), and g v ± (k) denote the subgradients of 
the dual function with respect to //(fc), A'(fc), and v f {k), 
respectively. These subgradients can be expressed in the fol- 
lowing simple forms 



g fit (k)=SR t -J2( P GZ- p Gj k )) 



(15a) 



A' 



g xt (k) = L* + Pk (*) + E P k W 



9=1 



A I 



E P k( k ) - P k(k) 



1 Ak)^P t R {k) 



J 

E p ^(fc)-^AW 

i=i 



(15b) 



(15c) 



where P^Jk), P*(fc), P^(fc), P* B .(k), P*(fc), and (fc) 
are given by (fT~6T> — (f20b- 

Iterations are initialized with arbitrary A*(0),z/(0) G M, 
and //(0) > 0. The iterates are guaranteed to converge to a 
neighborhood of the optimal multipliers l24l Sec. 6.3]. The 
size of the neighborhood is proportional to the stepsize, and 
can therefore be controlled by the stepsize. 

When the primal objective is not strictly convex, a primal 
averaging procedure is necessary to obtain the optimal power 
schedules, which are then given by 



J=0 



-x(fc-l). (21) 



{P G Jk)}U£ argmin j £ (c^J + (//(fc) - A'(ft))i$ ra ) j 

s.t. I8bl-t8dl 

{ p k( k )}J=i e argmin { £ f A*(fc)P^ - t/k(^j) j 

s .t. d3 

6 argmin { f] (a*^, - ^(P^)) ) 
U=i V / J 

s.t. (8i}-(D 

{^4 (*)}£=! e argmin | £ (^P^ + Fj(Pj) U 

s.t. ^-g5J 

{PA(fc),PA(fc)}Li£ argmin ( £ f („*(fc) - A^fc))^) + G({Pl» - £ ^(^1 
{PA s.t. (EJ, I t=i V / t=i J 



(16) 



(17) 



(18) 



(19) 



(20) 
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The running averages can be recursively computed as in ( |2TT >. 
and are also guaranteed to converge to a neighborhood of the 
optimal solution ||25ll . Note that other convergence-guaranteed 
stepsize rules (e.g., diminishing stepsize) and primal averaging 
(weights proportional to the stepsize) can also be utilized; 
see 1 26 1 for detailed discussions. 

2) Distributed Implementation: The form of the subgradi- 
ent iterations easily lends itself to a distributed implementation 
utilizing the control and communication capabilities of a 
typical microgrid. 

Specifically, the MGEM maintains and updates the La- 
grange multipliers via (TT4b . The LCs of conventional gen- 
eration, dispatchable loads, storage devices, and RES solve 
subproblems (fT6]i-(l20l), respectively. These subproblems can 
be solved if the MGEM sends the current multiplier iterates 
H*(k), A*0)> and v\k) to the LCs. The LCs send back 
to the MGEM the quantities Y%=\ P G m ( k )> E^=i p h n ( k )> 
E 9 Q =i^,(fc), T,U P bM ^d Pm which are 

in turn used to form the subgradients according to (TT3T >. The 
distributed algorithm using dual decomposition is tabulated as 
Algorithm Q] and the interactive process of message passing 
is illustrated in Fig. [2] 



quadratic. Correspondingly, the first four subproblems (fTSTi- 
([T9l > are essentially linear programs (LPs) or quadratic pro- 
grams (QPs), which can be solved efficiently. Therefore, the 
main focus is on solving (l20b . 

The optimal solution of P R (k) in d20l > is easy to obtain as 



nnun 
pmax 



if v\k) > A'(fc) 
if v\k) < A*(fc). 



(22) 



However, due to the absolute value operator and the maximiza- 
tion over w in the definition of G({P R }), subproblem ([20l > 
is a convex nondifferentiable problem in {P R }, which can be 
challenging to solve. As a state-of-the-art technique for convex 
nondifferentiable optimization problems l27l . 11241 Ch. 6], the 
bundle method is employed to obtain {P R (k)}. 
Upon defining 

T 

G({P R }) := G({P R }) - ]T ^(k)P R (23) 
t=i 

the subgradient of G({P R }) with respect to P R needed for 
the bundle method can be obtained by the generalization of 
Danskin's Theorem E4l Sec. 6.3] as 



Algorithm 1 Distributed energy management algorithm 



a 



Initialize Lagrange multipliers A* = p l = v l = 
repeat (k = 0, 1, 2, . . .) 
for t = 1,2,..., T do 

Broadcast A*(fc), and z/(fc) to LCs of con- 

vectional generators, controllable loads, storage devices, 
and RES facilities 

Update power scheduling Pq (k), Pjj (k), 
P| a (fc), P^.(fe), P R {k), and P R {k) by solving $&§2^ 

Update A*(fc), (/(k), and v\k) via di 
end for 

Running averages of primal variables via ( f2TT > 
until Convergence 




'pt pt 

R> R 



Fig. 2. Decomposition and message exchange. 

B. Solving the LC Subproblems 

This subsection shows how to solve each subproblem (fT~6b — 
©. Specifically, C^(-), -K n (-), -U* Eq (-), and Hj(-) are 
chosen either convex piece-wise linear or smooth convex 



A* — v l , if P R >J2(W{)* 



dG({P R }) = { (24) 

p -a* -i/*, if P t R <J2(wD* 



where for given {P R } it holds that 

w* G argmax ( ]T ^5*1^ - ^ Wf\ + 7*^ - E W*)H. 
wew I t=i V i=i i=i / J 

(25) 

With p := [P R ,...,P R ], the bundle method generates a 
sequence {p^} with guaranteed convergence to the optimal 
{P R (k)}; see e.g., E7), EU Ch. 6]. The iterate p e+1 is 
obtained by minimizing a polyhedral approximation of G(p) 
with a quadratic proximal regularization as follows 



pe+i := argmm 
pe 



infG,(p) 



(26) 



where G f (p) := max{G(p ) + go(p^- Po), • ■ • , G(p^) + 
gf(P — P<")}' §f i s th e subgradient of G(p) evaluated at the 
point p = pf, which is calculated according to ( f24b ; proximity 
weight is to control stability of the iterates; and the proximal 
center yg is updated according to a query for descent 



ye+i 



ye, 



if G(yt) - G(pt+i) >9ru 
otherwise 



(27) 



G(ye) - [G e (pe +1 



f\ 



Pe+i - ye\ 



where % 
(0,1). 

It is worth mentioning that d26| i is essentially a QP over 
a simplex in the dual space, which is efficiently solvable by 
practical optimization algorithms. The corresponding transfor- 
mation is shown in Appendix J] for the interested readers. 

Algorithms for solving (|25T l depend on the form of the 
uncertainty set W, and are elaborated next. 
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C. Vertex Enumerating Algorithms 

In order to obtain w*, the convex nondifferentiable function 
in d25l l should be maximized over W. This is generally an NP- 
hard convex maximization problem. However, for the specific 
problem here, the special structure of the problem can be 
utilized to obtain a computationally efficient approach. 

Specifically, the global solution is attained at the extreme 
points of the poly tope l24l Sec. 2.4]. Therefore, the objective 
in d25l l can be evaluated at all vertices of W to obtain 
the global solution. Since there are only finitely many ver- 
tices, (l25T l can be solved in a finite number of steps. 

For the poly topes W with special structure [cf. ©, ©J, 
characterizations of vertices are established in Propositions [3] 
and |4] Capitalizing on these propositions, vertex enumerating 
procedures are designed consequently, and are tabulated as 
Algorithms [2] and [3] 

Proposition 3. For a polytope A := {a G K"|a < a -< 
a, a mln < l'a < a max }, a v £ A is a vertex (extreme point) 
of A if and only if it has one of the following forms: i) a\ = 
flj or Hi for i = 1, . . . , n; or ii) a\ = a min —J2j i a a j or ftmax ~ 
Y,&i a j' a ) = % or f or i,j = l,---,n,j j^l 

Proof: See Appendix III-AI ■ 
Essentially, Proposition [3] verifies the geometric character- 
ization of vertices. Since W is the part of a hyperrectangle 
(orthotope) between two parallel hyperplanes, its vertices can 
only either be the hyperrectangle's vertices which are not cut 
away, or, the vertices of the intersections of the hyperrectangle 
and the hyperplanes, which must appear in some edges of the 
hyperrectangle. 

Next, the vertex characterization of a polytope in a Cartesian 
product formed by many lower-dimensional polytopes like A 
is established, which is needed for the uncertainty set @. 

Proposition 4. Assume b G R™ is divided into S consecutive 
and non-overlapping blocks as b = [b' 1; . . . , b' s ]', where b s G 
R' ls and E s =l n s — n - Consider a polytope B := {b G 
R n \b r< b r< b, bf in < V n b s < bf ax , s = 1, . . . , S}. Then 
b v = [(b^)', . . . , (bg)']' is a vertex of B if and only if for 
s = 1, . . . , S, h v s is the vertex of a lower- dimensional polytope 
B s := {b s G R n * \b s ±b s < b s , bf in < l' n b s < b™ x }- 

Proof: See Appendix IH-BI ■ 



Algorithm 2 Enumerate all the vertices of a polytope A 
l: Initialize vertex set V = 

2: Generate set A := {a G R"|5, = or a,, i = 1, . . n); 
check the feasibility of all the points in set A, i.e., if 

a min < < a roax| ; then V = V U {a} 

3: Generate set A := {a G M"|a ; = a min - 
Ej^aj or a max - E^iOj) & o = % or hO = 
l,...,n,j 7^ i}\ check the feasibility of all the points 
in set A, i.e., if a ^ a X a, then V = V U {a} 



Algorithms |2] and [3] can be used to to generate the vertices 
of uncertainty sets and (|6]l as described next, 
i) For uncertainty set (O, first use Algorithm [2] to obtain 
the vertices corresponding to each sub-horizon 7i jS for 



Algorithm 3 Enumerate all the vertices of a polytope B 
l: for i = 1,2, ... ,5 do 

2: Obtain vertex set V s by applying Algorithm |2] to B s 
3: end for 

4: Generate vertices b v for B by concatenating all the indi- 
vidual vertices b s as b v = [(hi)', (b v s )'}', b s G V s 



all the RES facilities. Then, concatenate the obtained 
vertices to get the ones for each RES facility by Step 4 
in Algorithm [3] Finally, run this step again to form the 
vertices of dUi by concatenating the vertices of each Wi. 
ii) For uncertainty sets ©, use Algorithm [2] to obtain the 
vertices for each sub-horizon T s . Note that concatenating 
step in Algorithm [3] is not needed in this case because 
problem (l25l l is decomposable across sub-horizons T s , 
s = 1 , . . . , S, and can be independently solved accord- 
ingly. 

After the detailed description of vertex enumerating proce- 
dures for RES uncertainty sets, a discussion on the complexity 
of solving d25b follows. 

Remark 2. (Complexity of solving d25ll). Vertex enumeration 
incurs exponential complexity because the number of vertices 
can increase exponentially with the number of variables and 
constraints l28l Ch. 2]. However, if the cardinality of each 
sub-horizon T s is not very large (e.g., when 24 hours are parti- 
tioned into 4 sub-horizons each comprising 6 time slots), then 
the complexity is affordable. Most importantly, the vertices of 
W need only be listed once, before optimization. 

IV. Numerical Tests 

In this section, numerical results are presented to verify the 
performance of the robust and distributed energy scheduler. 
The considered microgrid consists of M = 3 conventional 
generators, = 6 class tC\ dispatchable loads, Q = 6 class 
K-2 dispatchable loads, J = 3 storage devices, and 1 = 2 
renewable energy facilities (wind farms). The time horizon 
spans T = 8 hours, corresponding to the interval 4pm- 12am. 
The generation costs C m (Pc m ) = a mPc + b m Pc m and the 
utilities of class JC\ elastic loads U u (PdJ) = c nP[) n +d n Prj n 
are set to be quadratic and time-invariant. Generator param- 
eters are given in Table U and SR* = lOkWh. The relevant 
parameters of two classes of dispatchable loads are depicted 
in Tables |Il] and [ill] (see also 0, @, and J26)). The utility 
of class IC2 loads is (Pg ) := 7r*P|; with weights 
= 4, 3.5, . . . , 1, 0.5 for t = 4pm, . . . , 11am and q G Q. 

Three batteries have capacity 5™ ax = 30kWh (similar to 
13). The remaining parameters are Pj? ta = -lOkWh, P^. ax = 
lOkWh, Bf in = 5kWh, and if s = rf> = 0.9, for all j G 
J. The penalty weights for t = 4pm, . . . , 11am are set as 
w\ = 0.05, 0.1, . . . , 0.35, 0.4 and w\ = 0.1, j = 2, 3, t G T. 
The joint uncertainty model with S — 1 is considered for W 
[cf. ©], where Wf m = 40kWh, and Wj™* = 360kWh. In 
order to obtain V£- and W\ listed in Table |IV] MISO day 
ahead wind forecast data ||29l are rescaled to the order of 
1 kWh to 40 kWh, which is a typical wind power generation 
for a microgrid [16]. 



TABLE I 

Generating capacities, ramping power, and cost coefficients. 
The units of a m and b m are $/(kWh) 2 and $/kWh, respectively. 



Unit 


pmm 


pmax 


-^m.up(down) 


Clm 


bm 


1 


10 


50 


40 


0.006 


0.5 


2 


8 


45 


40 


0.003 


0.25 


3 


15 


70 


40 


0.004 


0.3 



TABLE II 

Class /Ci dispatchable loads parameters. The units of c„ and 
d n ARE $/(kWh) 2 and $/kWh, respectively. 



Ki 


Load 1 


Load 2 


Load 3 


Load 4 


Load 5 


Load 6 


pmm 

D n 
pmax 


0.5 
10 


4 
16 


2 
15 


5.5 

20 


1 

27 


7 

32 


Cn 
dn 


0.002 
0.2 


0.0011 
0.11 


0.003 
0.3 


0.0024 
0.24 


0.0015 
0.15 


0.0037 
0.37 



Similarly, the fixed load L l in Table IVl is a rescaled version 
of the cleared load provided by MISO's daily report 1301 . For 
the transaction prices, two different cases are studied as given 
in Table IVl where {a*} in Case A are real-time prices of the 
Minnesota hub in MISO's daily report. To evaluate the effect 
of high transaction prices, {a*} in Case B is set as 20 times 
of that in Case A. For both cases, /3* = 0.9a*, which satisfies 
the convexity condition for (PI) given in Proposition [TJ 

First, convergence of the Lagrange multiplier {A*} corre- 
sponding to the balance equation ( |8ot is confirmed for Case A 
by Fig. [3] It can be seen that A* converges for all i G T within 
a couple of hundred iterations, which verifies the validity of the 
dual decomposition approach using the subgradient projection 
method. With the running averages, convergence of the other 
dual and primal variables was also verified (for Case B too), 
but is omitted for brevity. 

0.4 1 . 1 1 1 1 1 



0.35- 




Iteration 

Fig. 3. Convergence of {A*}. 

The optimal microgrid power schedules of two cases are 
shown in Figs. [4] and [5] The stairstep curves include P G := 
Em P G m - p h = £„ Ph n , ™d P E = E 9 P l Eq denoting the 
total conventional power generation, and total elastic demand 
for classes K,\ and K.2, respectively, which are the optimal 
solutions of (P2). Quantity W^oret denotes the total worst-case 



TABLE III 

Class /C2 dispatchable loads parameters 



IC 2 


Load 1 


Load 2 


Load 3 


Load 4 


pram 
^E q 














pmax 


1.2 


1.55 


1.3 


1.7 


gmax 


5 


5.5 


4 


8 


k 


6pm 


7pm 


6pm 


6pm 




12am 


1 lam 


12am 


12am 



TABLE IV 
Limits of forecasted wind power 



Slot 


1 


2 


3 


4 


5 


6 


7 


8 


w\ 


2.47 


2.27 


2.18 


1.97 


2.28 


2.66 


3.1 


3.38 


w\ 


24.7 


22.7 


21.8 


19.7 


22.8 


26.6 


31 


33.8 


Wi 


2.57 


1.88 


2.16 


1.56 


1.95 


3.07 


3.44 


3.11 


wt 


25.7 


18.8 


21.6 


15.6 


19.5 


30.7 


34.4 


31.1 



wind energy at slot t, which is the optimal solution of (f25T > 
with optimal P R . 

A common observation from Figs. [4] and E] is that the total 
conventional power generation Pq varies with the same trend 
across t as the fixed load demand L*, while the K,\ elastic 
load exhibits the opposite trend. Because the conventional 
generation and the power drawn from the main grid are 
limited, the optimal scheduling by solving (P2) dispatches less 
power for P l D when L* is large (from 6PM to 9PM), and vice 
versa. This behavior indeed reflects the load shifting ability 
of the proposed design for the microgrid energy management. 
Due to the starting time S q (cf. Table HHI). another fact for both 
cases is the zero power scheduling of P l E for the time period 
of 4PM to 6PM, and the decreasing trend of P| after 6PM. 
Further elaboration on the behavior of each P E is provided 
later. 

Furthermore, by comparing two cases in Figs. [4] and [5] it 
is interesting to illustrate the effect of the transaction prices. 
Remember that the difference between P R and W^orst i s 
the shortage power needed to purchase (if positive) or the 
surplus power to be sold (if negative), Figs. [4] shows that 
the microgrid always purchases energy from the main grid 
because P R is more than Wl orst . This is because for Case 
A, the purchase price a* is much lower than the marginal 
cost of the conventional generation (cf. Tables HI andlVb. The 
economic scheduling decision is thus to reduce conventional 
generation while purchasing more power to keep the supply- 
demand balance. For Case B, since a* is much higher than that 
in Case A, less power should be purchased which is reflected 
in the relatively small gap between P R and W^, orst across time 
slots. It can also be seen that P R is smaller than W^ orst from 
7PM to 8PM, meaning that selling activity happens and is 
encouraged by the highest selling price /3* in this slot across 
the entire time horizon. Moreover, selling activity results in 
the peak conventional generation from 7PM to 8PM. Fig. [6] 
compares the optimal costs for the two cases. It can be seen 
that the optimal costs of conventional generation and worst- 
case transaction of Case B are higher than those of Case A, 
which can be explained by the higher transaction prices and 
the resultant larger DG output for Case B. 
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TABLE V 

Fixed loads demand and transaction prices. The units of of 

AND /?* ARE 0/KWH. 



Slot 


1 


2 


3 


4 


5 


6 


7 


8 


L* 


57.8 


58.4 


64 


65.1 


61.5 


58.8 


55.5 


51 


(Case A) 


















of 


2.01 


2.2 


3.62 


6.6 


5.83 


3.99 


2.53 


2.34 


P 


1.81 


1.98 


3.26 


5.94 


5.25 


3.59 


2.28 


2.11 


(Case B) 


















of 


40.2 


44 


72.4 


132 


116.6 


79.8 


50.6 


46.8 


P 


36.18 


39.6 


65.16 


118.8 


104.94 


71.82 


45.54 


42.12 



120 



100 



80 - 



■g 60 



% 40 



20 - 



■P'g 

•Pr 
■P'd 

P'e 
V 

WL 



4PM 5PM 6PM 7PM 8PM 9PM 10PM 11PM 12AM 
Time slot 



Fig. 4. Optimal power schedules: Case A. 



The optimal power scheduling of class AC 2 elastic load is 
depicted in Fig. [7] for Case A. The decreasing trend for all the 
AC 2 loads is due to the decreasing weights {ir^} from S q to 
T q , which is established from the fast charging motivation for 
the PHEVs, for example. 

Fig. |8] depicts the optimal SoC of the battery storage devices 
for Case A. The coincidence of SoC and the battery capacity 
for the several last slots is the result of the postulated battery 
cost penalizing the deviation from capacity. The increasing 
SoC behavior can be explained by the increasing weights 
H>. 

Fig. [9] shows the effect of different selling prices on 
the optimal energy costs, where Case B is studied with fixed 
purchase prices {a*}. It can be clearly seen that the net cost 
decreases with the increase of the selling-to-purchase-price 
ratio /3*/a*. When this ratio increases, the microgrid has a 
higher margin for revenue from the transaction mechanism, 
which yields the reduced worst-case transaction cost. Finally, 
it is worth mentioning that if a similar numerical test was 
implemented for Case A, but no significant change emerged 
in either type of cost. The reason is that the low real-time 
prices {a*} and in Case A do not considerably affect 
the transaction amount for different selling prices . 

V. Conclusions and Future Work 

A distributed energy management approach was developed 
tailored for microgrids with high penetration of renewable 
energy sources. By introducing the notion of committed 
renewable energy, a novel model was introduced to deal 




4PM 5PM 6PM 7PM 8PM 9PM 10PM 11PM 12AM 
Time slot 



Fig. 5. Optimal power schedules: Case B. 
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I Conventional generation cost 
| Negative utility of K loads 

| Negative utility of K 2 loads 

] Battery cost 

I Worst-case transaction cost 
] Microgrid net cost 



o 
o 

| 100- ■ 
50- I 

^| 

-50 - 



Case A 



Fig. 6. Optimal costs: Case A and B. 



Case B 



with the challenging constraint of the supply-demand balance 
raised by the intermittent nature of renewable energy sources. 
Not only the conventional generation costs, utilities of the 
adjustable loads, and distributed storage costs were accounted 
for, but also the worst-case transaction cost was included in 
the objective. To schedule power in a distributed fashion, the 
dual decomposition method was utilized to decompose the 
original problem into smaller subproblems solved by the LCs 
of conventional generators, dispatchable loads, DS devices and 
the RES. The numerical tests demonstrated the optimal power 
generated, consumed by and delivered to the DERs, and the 
controllable loads across the entire time horizon. 

A number of interesting research directions open up towards 
extending the model and approach proposed in this paper. 
Some classical but fundamental problems, such as the optimal 
power flow (OPF) and the unit commitment (UC) problems 
are worth re-investigating with the envisaged growth of RES 
usage in microgrids. 
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Fig. 7. Optimal power schedule for Pi : Case A. 



35 
30 



I ft, 
\Pb, 

]Pb, 



s 25 

3 

oi 20- 
3 

15- 

o 
w 

™" 10- 
5 - 



II II 



4PM 5PM 6PM 7PM 8PM 9PM 1 0PM 1 1 PM 1 2AM 
Time slot 



Fig. 8. Optimal power schedule for Case A. 



Appendix I 
Enhancing the Bundle Method 



, Pi ll ll2 

mm r+— \\p - ye\\ 



Using an auxiliary variable r, d26l i can be re-written as 

(28a) 

s.t. G( Pi ) + g<(p - Pi ) < r, i = 0, 1,. . .,£. (28b) 

Introducing multipliers £ 6 K.^ 4 " 1 , the Lagrangian is given as 

t+i 



P l \\ l|2 



\ i=0 ' 

e+i 

+ X)Ci(G(p0+g{(p-Pi)). (29) 



i=0 



Optimality condition on p, i.e., V p £(r, p, £) = , yields 



P = ye 



1 <+l 

— y^&gi- 

pi 



(30) 
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Fig. 9. Optimal costs: Case B. 



Substituting OOJl into 



max 



e+i 



i=0 



the dual of 

e+i 



IS 



s.t. 



I h 0, 1'^ = 1 



^i{G{pi)+^(yt-Pi)) (31a) 

(31b) 



i=0 



where 1 is the all-ones vector. 

Note that OTb is essentially a QP over the simplex in K +1 , 
which can be solved very efficiently. 

Appendix II 
Proofs of Propositions 

To prove Propositions [3] and [4] the following lemma is 
needed, which shows sufficient and necessary conditions for 
a point to be a vertex of a polytope represented as a linear 
system (31] Sec. 3.5]. 

Lemma 1. For a polytope V := {x e R™|Ax < c}, a point 
v G V is a vertex if and only if there exists a subsystem 
Ax < c of Ax < c so that rank(A) = n and v is the unique 
(feasible) solution of Av = c. 



A. Proof of Proposition [3] 

The polytope A := {a G R"|a < a ^ a, a min < l'a < 
} can be re-written as A := {a G M™|Aa ^ c}, where 



1, —1]' and c := [a' 



a 
A 

By Lemma [T] enumerating vertices of A is equivalent to 
finding all feasible solutions of the linear subsystems Aa = c, 
such that rank-n matrix A is constructed by extracting rows 
of A. It can be seen that such full column-rank matrix A can 
only have two forms (with row permutation if necessary): i) 
Ai = diag(d) with di G { — 1, 1}, i = 1, • • • , n; ii) A.2{i, ■) = 
±1', i = 1, . . . , n, and A.2(j, :) = Ai(j, :), Vj ^ i. Basically, 
Ai is constructed by choosing n vectors as a basis of R" from 
the first 2n rows of A. Substituting any row of A x with ±1', 
forms A 2 . Finally, by solving all the linear subsystems of the 
form Afca = c^, for k = 1,2, Proposition [3] follows readily. 
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B. Proof of Proposition __ 

The polytope B := {b E M"|b ^ b ^ b, bf in < l' n b s < 
6™ ax , s = 1, . . . , S} can be re-written as 6:={be K n |Bb -< 
c}, where B := diag(B 1; . . . , B s ), c := [ci, . . . , c' s ]', B s := 

[Ira s xn„ j — In, xn s 7 1; — an d c s : = [ a s , — &s , &™ ax , — 6™ ln ] ' 

for s = 1,...,S. 

Similarly by Lemma __ all the vertices of B can be enumer- 
ated by solving Bb = c, where the rank-n matrix B is formed 
by extracting rows of B. Due to the block diagonal structure 
of B, it can be seen that the only way to find its n linear 
independent rows is to find n a linear independent vectors from 
the rows corresponding to B s for s = 1,...,S, In other 
words, the vertices b v can be obtained by concatenating all 
the individual vertices h s as stated in Proposition __ 
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